On 
On 
On 






-<— > 

CO 



Density Matrices for a Chain of Oscillators 



Ingo Peschel* and Ming-Chiang Chung^ 



Fachbereich Physik, Freie Universitat Berlin, 
Arnimallee 14, D-14195 Berlin, Germany 

February 1, 2008 



C3 . Abstract 

s. 

i I We consider chains with an optical phonon spectrum and study the reduced 

density matrices which occur in density-matrix renormalization group (DMRG) 
calculations. Both for one site and for half of the chain, these are found to be 
O | exponentials of bosonic operators. Their spectra, which are correspondingly expo- 

nential, are determined and discussed. The results for large systems are obtained 

C*~) . from the relation to a two-dimensional Gaussian model. 

> 

(N ! 1 Introduction 

\o 

The success of the density-matrix renormalization method (DMRG) in treating one- 
dimensional quantum systems [[l], Q is closely related to the properties of the involved 
density matrices. In the procedure, one determines the eigenvectors of these matrices and 
uses those with the largest eigenvalues as a truncated basis. To be able to single out a 
relatively small number, however, the density-matrix spectrum has to decrease rapidly 
enough. Indeed, it is usually found in the numerical calculations that the eigenvalues 
decay roughly exponentially. 

In a previous publication || it was pointed out that, for non-critical integrable models, 
the exponential behaviour is ultimately a consequence of the Yang-Baxter equations. For 
two spin one-half models, the transverse Ising chain and the uniaxial Heisenberg chain, 
analytical formulae were given and verified in detail in DMRG calculations. 

In the present article, we want to extend these considerations to phonons, i.e. to a 
bosonic problem. So far, comparatively few DMRG studies have been dealing with bosons 
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||]-[[nj. The difference to spin systems is that the full Hilbert space always has infinite 
dimension. Therefore any numerical treatment has to start with a truncation. One can 
do this in analogy with the DMRG procedure by selecting local states via the density 
matrix for a single site || |4|]. This is still a nontrivial quantity with an infinite number 
of eigenstates in a full treatment, and it is of interest to find its properties in a solvable 
case. The same holds, of course, for the more complicated density matrix of a half-chain 
which is used in the DMRG algorithm. 

The system we study is a purely bosonic model, a chain of L harmonic oscillators 
with frequency u , coupled together by springs. It has a gap in the phonon spectrum 
and is a non-critical integrable system as the spin models mentioned above. We write the 
Hamiltonian 

H = E(-^+^o^) + E^(^ + i-^) 2 (i) 

and will frequently use the form u = 1 — k, so that for k = there is no dispersion, 
while for k — > 1 the spectrum becomes acoustic and the system critical. 

We first consider in Section |2] the density matrix p\ for one oscillator and show that 
it can always be written as the exponential of the Hamiltonian of a (new) harmonic 
oscillator. The spectrum therefore is purely exponential, with a decay rate depending 
on k and (weakly) on the chosen site. This generalizes the known result for the case 
L = 2 |14|| . The eigenfunctions have the character of squeezed states and are used later 
for numerical calculations. In Section [|, we turn to the density matrix ph for half of the 
system. We treat the case of small and large L explicitely and find that ph has the same 
exponential form, with the number of oscillators in the exponent determined by the size 
of the system. The result in the thermodynamic limit is derived by relating the chain 
to a massive two-dimensional Gaussian model and its corner transfer matrices (CTMs). 
It is very similar to that for the spin chains in || which lead to fermionic operators 
instead of bosonic ones. In particular, the spectrum without the degeneracies is purely 
exponential. Its form for different values of k and different sizes L is discussed in more 
detail in Section f|, including numerical results obtained by truncation and by DMRG 
calculations. These also illustrate, to what extent the degeneracies are reproduced in an 
approximate treatment. The concluding Section [| contains a summary and additional 
remarks. Some details concerning the case L = 4 and the Gaussian model are given in 
two appendices. 

2 Density Matrix for One Oscillator 

In this section we consider the case where one oscillator is singled out and all other 
form the environment. The corresponding density matrix (determined numerically) was 
used previously in the study of an electron-phonon system M. Here, it can be found 



analytically. 

The ground state of H in ([!]) has the form 

tt(sc) = C ■ exp(--J2 AjXiXj) (2) 

where x = (x\,X2, ■ • • , xi). The matrix 

A ij = £ u Mi)4>q{3) (3) 

is determined by the frequencies u; g and the eigenvectors <p q (i) of the normal modes. From 
the total density matrix 

p(x,x') = #(#)*(« (4) 

one then obtains the reduced one for oscillator I by integrating over all other coordinates 
Xi = x\. This leads to 

pi(xi,x[) = d ■ exp(--(a-6)x 2 ) exp(--(xi - x[) 2 ) exp (--(a - 6)xf ) (5) 



with the constants 



a = A u (6) 

6 = £ A u i^]- 1 A 3l (7) 



*jV 



where A^ is the matrix obtained from A by deleting the Z-th row and column. The second 
exponential in ([5]) can be transformed into a differential operator, giving 

1 1 d 2 1 

px = C 2 • exp (--^V ) exp (- — ) exp (--wV ) (8) 

where y 2 = bxf/2 and a; 2 /4 = (a — &)/&. Writing this in terms of Bose operators a , 
a^ one can bring it into diagonal form by an equation-of-motion method. The necessary 
Bogoljubov transformation is 

P = cosh 9- a + sinh0-cr (9) 

with 

e e = (1 + ^) 1/4 (10) 

As a result, one finds that px has the form 

px = K ■ exp (-ft) (11) 




Figure 1: Eigenvalue e in the density matrix for an oscillator at the end of a chain, as a 
function of k for different lengths L and Uq = 1 — k 



where 

H = eft (3 

is the Hamiltonian of a harmonic oscillator with energy 



2sinh 



-i 



.Co> 



2sinh~ 1 [\la/b 



1 



(12) 
(13) 



Therefore the eigenvalues of p\ are w n = Ke~ £n ,n > 0, and the spectrum is purely 
exponential. The constant K follows from the sum rule Tr(pi) = J2n w n = 1- 

This result is completely general. The details of the oscillating system and the position 
of the chosen oscillator enter only via the ratio a/b. The same constants a and b also 
determine the probability to find a certain elongation x\. However, as seen from pi(xi, xi) 
in (H), this quantity depends on the difference (a — b) and thus has no direct relation to 
e. 

In the simplest case of two oscillators (L = 2) one finds explicitely 



or, equivalently, 



e = 2 sinh W AuiUzj '(u>i — a; 2 ) 2 



ln(coth 2 (^)) 



(14) 



(15) 
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Figure 2: Eigenvalue e as a function of the position of the oscillator, for a chain of L = 30 
sites and three different values of k. 



Uq + 2k are the two eigenfrequencies and e 2v = uj 2 /uji- This is 



yj in a different way. 



where U\ = uj , uj 2 = 
the result obtained in 

In Fig.^, e is shown as a function of k, putting uj = (1 — k). For k — > it diverges 
logarithmically. In this limit the influence of the second oscillator vanishes, *£?(x) becomes 
a product state and one is left with only one nonzero eigenvalue Wq — 1. For k — > 1, on 
the other hand, e goes to zero as \/I — A; and the eigenvalues w n decrease only very slowly, 
which reflects the strong coupling. These features are encountered also in all other cases. 
For L = 3 one can still give explicit analytical formulae, but for larger L the problem 
has to be treated numerically. In the figure, two additional cases, L = 10 and L = 100, 
are shown. The limit L — > oo, which is approached exponentially in L with a correlation 
length increasing with k, is indistinguishable from L = 100 on the given scale. 

One can also investigate how s varies with the position along the chain. The result 
is shown in FigfJ for several values of k. One sees that e is large at the ends. This 
corresponds to the fact that the influence of the environment is smaller there. At the 
next site, however, e drops and then approaches the bulk value from below as one moves 
into the interior. The approach becomes slower as k increases. The overall differences in 
the e-values are not very large, though, as seen in the figure. 

Due to the form of p±, its eigenstates are standard oscillator functions of a coordinate z 



which differs from x\ by a scale factor. Compared with the eigenfunctions of the uncoupled 
oscillator I, they are squeezed states whose spatial extent is reduced by a factor q = Wcuo/7, 

where 7 : 



'a(a — b). For small k, q approaches one and the two sets of functions coincide. 
With increasing k, the amount of squeezing increases, and it is then advantageous to 
choose the squeezed states as a local basis. This was done in the numerical calculations 
which will be presented in Section [|. 



3 Density Matrix for a Half-Chain 

We now turn to the central quantity in the usual DMRG calculations, the reduced density 
matrix for half of the system. It enters each time the system is enlarged in the infinite-size 
algorithm. We will determine its spectrum in the two limits of small and large L. 

For L = 2, one-half of the system is just one oscillator and ph has already been 
obtained in Section H We therefore proceed immediately to the case L — 4. First, we 

1/2 

note that the square root p h follows directly from \I/. If the coordinates along the chain 
are (x 2 ,Xi,x' 1 , x' 2 ), one has 



1/2 



Ph {Xi,X2',X l ,X 2 ) — ty(X2,Xi,X 1 ,X 2 ) 

Taking into account the form (0) and the symmetries, this leads to 



(16) 



1/2 
Ph 



C ■ exp < -- Y^ aij(xiXj + x'ix'j) - Y b 



ijXiXj 



(17) 



1.1 



where the symmetric (2 x 2) matrices a^ and b^ follow from the matrix A of Section |^. 
Altogether one has six different coefficients which couple the variables as shown in the 
following figure. 
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Figure 3: Lowest eigenvalues in the density matrix of a half-chain. Plotted are e for 
L = 2 and E\ for L = 4 and L = oo. 



The cross-couplings, shown as dashed lines, can be eliminated by introducing new 
coodinates y^ y[. After that, a sequence of transformations similar to those in Section |] 
brings p h (and thus ph) into diagonal form. Some of the details are given in Appendix 0. 
The final result is that ph has also the form (O) where now 



n = E £$& 



(« 



i=i 



describes two harmonic oscillators with energies ej and £2- Thus one obtains a simple 
generalization of the case L = 2. Also the variation of £\ with fc is very similar to that of 
e in Section ||. This is shown in Fig||, where both quantities are plotted. In particular, 
one finds that they coincide in the limit k — > 0. The ratio £2/^1 equals 3 for small k, 
drops to a minimum of 2.866 for k = 0.34 and then increases continuously, because £2, in 
contrast to S\, stays finite as k — > 1. The shape of the spectrum, which depends on the 
ratio £2/^1, w iH be discussed in the next section. 

At this point one can already conjecture that the structure of ph remains the same also 
for larger L. A direct computation as above does not seen feasible, though. In the limit 
of large L, however, a different approach is possible. As in |L7, ||| one first relates ph to 
the partition function of a two-dimensional classical system, which is a massive Gaussian 



model in our case, in the form of an infinite strip of width L with a perpendicular cut. 
This connection is discussed in more detail in Appendix |B[ One then expresses the 
partition function as the product of four corner transfer matrices. In the case where L 
is much larger than the correlation length, one can use the thermodynamic limit of these 
CTMs and finds for ph the form ( [TT]) with an operator 7i which is very similar to H in 
(H). The coefficients, however, are multiplied by additional site-dependent factors which 
increase linearly along the chain and reflect the corner geometry. Up to a prefactor it 
is the operator given in fl32|) in Appendix |B| and its diagonalization amounts to finding 
the normal modes of the corresponding vibrational problem. From the results in |1| one 
obtains 

n = E(2j-i)^]A- (19) 

with 



where I(k) is the complete elliptic integral of the first kind and k' = yl — k 2 . Therefore 
Ti describes an infinite set of harmonic oscillators with energies Ej = (2j — l)e and is a 
straightforward extension of the results for small L. 

The parameter e = e\ is also shown in Fig|J. For k — > 0, it has exactly the same 
expansion as e for L = 2 and e± for L = 4. For k — > 1, it vanishes only logarithmically, 
i.e. more slowly than the quantities for finite L. 

One should note that the results (|19|) , (p0|) are formally the same as for the transverse 
Ising chain in the disordered phase ||. The only difference is that there the operators 
/3,/3^ are fermionic (so that ft (3 = 0, 1), whereas here they are bosonic. Such similarities 
can also be observed in the row transfer matrices of the Gaussian and the Ising model, if 



one uses the corresponding parametrizations [19|]. The consequences for the spectrum of 
Ph are discussed below. 

4 Spectra and Numerics 

In the following, we show the density-matrix spectra for half-chains and discuss some 
numerical aspects. In the figures, the eigenvalues w n of ph are ordered according to 
magnitude and plotted on a semilogarithmic scale. 

Figure [| shows spectra for L = 4 and several values of k. These results were ob- 
tained by calculating the two energies £i, £2 numerically from the formulae in Appendix JS]. 
Apart from the rapid decrease, one notes a clear ladder structure for the smallest three 
fc's. It results from the relation 62 — 3ei which leads to the approximate degeneracies 
(1, 1, 1, 2, 2, 2, 3) for the first seven levels. The steps for k = 0.3 are less perfect, since €2 
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Figure 4: Density-matrix spectrum for L = 4 and five different values of k. 



deviates more from 3ei in this case. For the two largest fc's, £2 — 4ei and £2 — 6ei, so 
that the first step appears at these levels and the spectra look more stretched out. 

It is interesting to see, how these results are recovered in a numerical treatment using 
a truncated Hilbert space. If one works with the eigenstates of pi, a small number (5 — 7) 
is sufficient for not too large k. For example, if k = 0.5 and one chooses the same r states 
(with some average e-value) for all four sites, the error in the ground-state energy E Q /L 
is of the order 10 _r . The spectra which then result are shown in Fig.|5] for three values of 
r. The first w n are always quite accurate, but there are characteristic differences for the 
following ones, which are connected with the number of steps, i.e. with the degeneracies. 
One can see that if r states are kept, the pattern is correct for the first r levels (counted 
from the top). At the next level, the state with energy r£\ is missing and the corresponding 
step is absent. Thus there is a certain correspondence between the states in the local basis 
and in the density matrix. For smaller w n , however, the situation is less clear, and the 
spectrum finally becomes irregular. The tails of the approximate spectra always lie below 
the exact one. 

In order to obtain results also for L > 4, we have carried out DMRG calculations, 
using 7 states at each site, with an e corresponding to L — 30. With m = 7 kept states 
per block, the error in Eq/L was about 3 ■ 10~ 7 for k = 0.5. Fig.|6] shows the resulting 
spectra for L = 6 and L = 14, together with the thermodynamic limit according to 
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Figure 5: Density-matrix spectrum for L = 4 and k = 0.5, calculated with different 
numbers of oscillator states. 



(0),(^3). One notes, that the spectra for the two L's are similar, though not identical. 
Compared to L = 4, the degeneracies have changed to (1, 1, 1,2,2,3,4). The latter two 
result from a third energy e 3 = 5ei, which first appears for L = 6. This shows that, 
indeed, the number of oscillators in p^ is equal to the size L/2 of the half-chain. One also 
sees that for L = 14 the first two steps have become perfect, so that e 2 = 3s i as for the 
infinite system. Up to some small deviations, this also holds for the next two steps. Only 
for the remaining levels 8 and 9, the degeneracies are not correct. This is the same effect 
as found above for L = 4. 

For L = 14, the Sj are also numerically very close to the large-L limit. For example, 
£i agrees with the exact result 4.0189 up to three decimal places. This can be understood 
from the short correlation length £ ~ 3 for k = 0.5 which makes size-effects small. Finally, 
we want to mention that, in the thermodynamic limit, the multiplicities are just one-half 
of those found in the fermionic case for the ordered phase where Sj = 2ej 0. This is 
because the number Pj of partitions without repetition is the same as that of the odd 
integers with repetition, Pj = PL-v Therefore the degeneracies for the bosonic case are 
not as large as one might expect at first. 
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Figure 6: Density-matrix spectrum for k = 0.5 and two sizes L, calculated with DMRG 
using 7 states and m = 7. Also shown is the analytical result for L — > oo. 



5 Conclusion 

We have investigated a bosonic system, where the ground-state density matrices can 
be determined explicitely in various cases. It turned out that they are exponentials 
of oscillator Hamiltonians, so that all results are quite transparent. The spectra have 
exponential character and the eigenf unctions are oscillator states. For the single-site 
density matrix, these states are related to those of the chain oscillators by squeezing. For 
the half-chain density matrix, they are connected with certain normal modes concentrated 
near the middle of the system. The thermodynamic limit was obtained in the same way 
as for the integrable spin chains treated previously, and the spectra are very similar to 
those found there. By counting the degeneracies, one would arrive at formulae as given 
in 
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Taking all this together, the chain treated here may serve as a standard example where 
one can see the features of the density matrices in detail. In this context, it would still be 
interesting to determine the half-chain density matrix for arbitrary sizes, in particular at 
the critical point, where the vibrational spectrum becomes acoustic. This case has already 
been studied by DMRG [|| but, as for the critical spin models, the density-matrix spectra 



11 



have yet to be explained. Another question is, whether the model of coupled oscillators, 
for which the ground state is known explicitely, could be used to study density matrices 
in higher dimensions. For the DMRG method, it would be quite important to know if the 
spectral properties change in this case. 
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Appendices 

A Four Oscillators 

In order to diagonalize p h in ([17|) one proceeds as follows. First, new coordinates are 
introduced by a rotation (xi,x 2 ) — > (2/1,2/2) with angle <p and analogously for the primed 
quantities. This leads to new quadratic forms in the exponent, with coefficients a^ and bij. 
Choosing tan2y? = 26 12 /(6 n — 622), the cross-term & 12 becomes zero. One then considers 
the factors 

exp (-1/2 a ti y\ ) exp (-b u y^ ) exp (-1/2 a u yf ) , i = 1, 2 (21) 

which contain only {yi-,y'j)- These can be transformed as in Section |2] and one obtains 
exponentials of harmonic oscillators with energies 

Vi = 2sinh- 1 (fi i /2) (22) 

where 



n t /2 = V / («n + ^)/(-2 6 ii ) (23) 

In terms of the new coordinates z* one then has 



Ph 2 = C ■ exp (-/1Z1Z2) exp (- ^(--7^2 + tAzI ) ) exp (-pz lZ2 ) (24) 



Here Zi = yi/\, \i = O12 AiA 2 and the A; are given by 



A - is 1+ f (25) 



'«"! , 



In the final step, one expresses (|2"4]) in terms of bosonic operators at, a} and considers 
Heisenberg-like operators p h ctip h which are found to be linear combinations of the 
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4- Therefore, a transformation as in the analogous fermionic case |16| 

Pi = J2(9ji&i + hjia\) 



(26) 



brings p h into the form (|TT|) , ([18|) with energies £j/2. These energies follow from a simple 
quadratic equation, namely 



cosh — 
2 



1 1 

-(ci + c 2 ) ±\/t( c i - c 2 ) 2 + 4p 2 sis 2 



(27) 



where q = cosh z/j, S; = sinh z/j andp = \xj2^/v\V2- 

These quantities have to be evaluated, starting from the initial constants o^ and bij, 
which are simple analytic expressions involving the four eigenfrequencies of the chain. It 
turns out that, for small values of k, the p-term in (B7I) is unimportant, which leads to 
£j = 2i/j and e 2 — 3ei. The ratio £2/^1 thus has the same value as in the thermodynamic 
limit. Moreover, E\ has the same asymptotic form, t\ = 21n(4//c), as for L = 2 and 
L — > 00. This can be attributed to the short correlation length which suppresses size 
effects in this limit. 



B Relation to the Gaussian Model 

The Hamiltonian H in ([!]) has a close relation to the transfer matrix of a two-dimensional 
Gaussian model (GM). The connection is the same as between the transverse Ising chain 
and the two-dimensional Ising model 0. Consider a lattice with variables x (—00 < x < 00) 
at each site, a nearest-neighbour coupling energy | K (x—x 1 ) 2 and an on-site energy | A x 2 , 
all in units of fcgT '. If the lattice is oriented diagonally, the appropriate transfer matrix 
T involves the piece shown in the figure below. 




Xi-i 



* 1 



Xi+l 



One can then verify by a simple direct calculation (using two interpenetrating lattices) 
that, with periodic boundary conditions, 



[H,T] 







(28 
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provided that k = K 2 and Uq = A(A + 4K). In this case, T and H have common 
eigenfunctions and \1/ in (H) gives the maximal eigenvalue for T. This allows one to obtain 
\1/ and also ph from the partition function of a two-dimensional system B I7|. If the GM 



has open boundaries, one has to modify H at the end, so as to preserve fl28|) . However, 
for a system with L ^> £, where £ is the correlation length given by £ = 2/ In (1/fc), this 
effect is not important and can be neglected. 

An alternative approach is to consider a Gaussian model with anisotropic couplings 
for periodic boundary conditions, to show that the T's for different anisotropics commute 
and to realize that a proper derivative leads to H |20|| . To do this, one uses an elliptic 
parametrization, so that the two couplings are, for example, 

K\ = —i/sn(iu,k) ; K 2 = iksn(iu,k) (29) 

with the Jacobi function sn of modulus k. This parameter also determines the on-site 
energy A and thus the distance to the critical point A = 0, as well as the correlation 
length. The parameter u, on the other hand, specifies the ratio K1/K2. It varies between 
and I{k'), where I is the complete elliptic integral of the first kind and k! = \/l — k 2 . 
The isotropic case corresponds to u = I(k')/2. (Our notation differs slightly from that 
in [20]. We have interchanged k <->■ k', written u instead of \9, used x = VA0 for the 
Gaussian variables and have set a = — 1.) The derivative (dlnT/du) then leads again 
to (|l]) with uiq = (1 — k), which is the reason for choosing this parametrization in H. 

As discussed in , the density matrix p h for half of the system is, for L ^> £ and up 
to a prefactor, 

p h = ABCD (30) 

where A, B, C, D are the corner transfer matrices for the four infinite quadrants of the two- 
dimensional system. Due to the integrability of the Gaussian model, i.e. the Yang-Baxter 
equations, these have the exponential form 

A = e ~ uHc ™ (31) 

and similarly for B,C,D, with TCctm given by 

Hctm = E {-^( 2 ^ - !)^ + l( 2n -!)(!- k ? x l + fnk(x n+1 - x n ) 2 \ (32) 

n nnnnppfiATi tsrifn f ne T-TamiUnni an hmif n — i. nf A 



This operator was studied in [JTHj] in connection with the Hamiltonian limit u — > of A 
where one can determine its form simply by inspection. It is associated with a corner of 
Ramond type, i.e. without a site at the tip. In terms of vibrations, it describes a system 
of coupled oscillators, where the spring constants and inverse masses increase along the 
chain. It can be diagonalized with the help of Carlitz polynomials and then becomes the 
sum of harmonic oscillators with eigenvalues (2j — l)ir/2I(k'). For ph one needs ABCD, 
or A A if one has an isotropic model. In either case this gives a factor 21 (k'), so that the 
energy becomes Ej = (2j — l)nl(k f )/l(k) and one arrives at the result fllTf) , (|20|) . 

14 



References 



[i] 

[2] 

[3] 
[4] 
[5] 

[6] 

[7 

[8 

[9 

[10 

[11 
[12 

[13 

[14 

[15 



[16 

[17; 



S.R. White, Phys. Rev. Lett. 69, 2863 (1992); Phys. Rev. B 48, 10345 (1993) 

for a review see: Density-Matrix Renormalization, I. Peschel, X. Wang, M. Kaulke 
and K. Hallberg(eds.), Lecture Notes in Physics Vol. 528, Springer (1999) 

I. Peschel, M. Kaulke and O. Legeza, Ann. Physik (Leipzig) 8, 153 (1999) 

for a brief review, see E. Jeckelmann, C Zhang and S.R. White in Ref. || 

L.G. Caron, S. Moukouri, Phys. Rev. Lett. 76, 4050; Phys. Rev. B 56, R8471 (1997) 

R.V. Pai, R. Pandit, H.R. Krishnamurthy and S. Ramasesha, Phys. Rev. Lett. 76, 
2937 (1996) 

E. Jeckelmann, S.R. White, Phys. Rev. B 57, 6376 (1998) 

C. Zhang, E. Jeckelmann and S.R. White, Phys. Rev. Lett. 80 , 2661 (1998) 

T. Kuhner, H. Monien, Phys. Rev. B 58, R14741 (1998) 

R.J. Bursill, R.H. McKenzie and C.J. Hamer, Phys. Rev. Lett. 80, 5607 (1998); 
|cond-mat/98T2409| (1998) 



R.J. Bursill, |cond-mat/98l2349| (1998) 



S. Rapsch, U. Schollwock and W. Zwerger, Europhys. Lett. 46 5 (1999) 

T.D. Kuhner, S.R. White and H. Monien, |cond-mat/9906019| (1999) 

D.Han, Y.S. Kim and M.E. Noz, Am. J. Phys. 67, 61 (1999) 

K. Okunishi, Y. Hieida and Y. Akutsu, Phys. Rev. E 59, R6227 (1999); |cond 



mat/98l0239| (1998) 



E.H. Lieb, T.D. Schultz and D.C. Mattis, Ann. Phys. (NY), 16, 406 (1961) 

T. Nishino, K. Okunishi, Density Matrix and Renormalization for Classical Lattice 
Models in: Strongly Correlated Magnetic and Superconducting Systems, ed. G.Sierra 
and M.A.Martin-Delgado, Lecture Notes in Physics Vol. 478, Springer Berlin, Hei- 
delberg (1997) (see also |cond-mat/9610107|) . 



[18] I. Peschel, T.T. Truong, Ann. Physik (Leipzig) 48, 185 (1991) 



15 



[19] M. Sato, J. Miwa and M. Jimbo : Holonomic Quantum Fields V, Publ. Res. Inst. 
Math. Sci. (Kyoto Univ.) 16, 531 (1980). 

[20] G.M. Babudjan, M.G. Tetelman, Theor. Math. Phys. 51, 484 (1982) 



16 



